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Abstract 

The interest of compressive sampling in ultrasound imaging has been recently extensively evaluated by several 
research teams. Following the different application setups, it has been shown that the RF data may be reconstructed 
from a small number of measurements and/or using a reduced number of ultrasound pulse emissions. Nevertheless, 
RF image spatial resolution, contrast and signal to noise ratio are affected by the limited bandwidth of the imaging 
transducer and the physical phenomenon related to US wave propagation. To overcome these limitations, several 
deconvolution-based image processing techniques have been proposed to enhance the ultrasound images. In this 
paper, we propose a novel framework, named compressive deconvolution, that reconstructs enhanced RF images 
from compressed measurements. Exploiting an unified formulation of the direct acquisition model, combining 
random projections and 2D convolution with a spatially invariant point spread function, the benefit of our approach 
is the joint data volume reduction and image quality improvement. The proposed optimization method, based on 
the Alternating Direction Method of Multipliers, is evaluated on both simulated and in vivo data. 

Index Terms 

Compressive sampling, deconvolution, ultrasound imaging, alternating direction method of multipliers 


1. Introduction 

U ltrasound (US) medical imaging has the advantages of being noninvasive, harmless, cost- 
effective and portable over other imaging modalities such as X-ray Computed Tomography or 
Magnetic Resonance Imaging [[1]. 

Despite its intrinsic rapidity of acquisition, several US applications such as cardiac, Doppler, elastog- 
raphy or 3D imaging may require higher frame rates than those provided by conventional acquisition 
schemes (e.g. ultrafast imaging [^) or may suffer from the high amount of acquired data. In this context, 
a few research teams have recently evaluated the application of compressive sampling (CS) to 2D and 3D 
US imaging (e.g. [|3]-[8l) or to duplex Doppler Q. CS is a mathematical framework allowing to recover, 
via non linear optimization routines, an image from few linear measurements (below the limit standardly 
imposed by the Shannon-Nyquist theorem) [fTOl ITT]| . The CS acquisition model is given by 

y = + n (1) 

where y G corresponds to the M compressed measurements of the image r G (one US 
radiofrequency (RF) image in our case), $ G represents the CS acquisition matrix composed for 

example of M random Gaussian vectors with M « N and n G stands for a zero-mean additive 
white Gaussian noise. 

The CS theory demonstrates that the N pixels of image r may be recovered from the M measurements 
in y provided two conditions: i) the image must have a sparse representation in a known basis or frame 
and ii) the measurement matrix and spasifying basis must be incoherent W2\ . In US imaging, despite the 
difficulties of sparsifying the data because of the speckle noise, it has been shown that RF images may 
be recovered in basis such as 2D Fourier 0|, wavelets waveatoms [[71 or learning dictionaries [IT3l . 
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considering Bernoulli Gaussian [IT4l or a-stable statistieal assumptions [|T5l and using various aequisition 
sehemes such as plane-wave IH, Xampling [[5l or projections on Gaussian ^ or Bernoulli random veetors 

0 . 

However, the existing methods of CS in US have been shown to be able to recover images with a 
quality at most equivalent to those acquired using standard sehemes. Nevertheless, the spatial resolution, 
the signal-to-noise ratio and the eontrast of standard US images (r in eq. ([T])) are affeeted by the limited 
bandwidth of the imaging transdueer, the physieal phenomenon related to US wave propagation sueh as 
the diffraetion and the imaging system. In order to overeome these issues, one of the researeh traeks 
extensively explored in the literature is the deeonvolution of US images [fT6] - l2T]l . Based on the first 
order Born approximation, these methods assume that the US RF images follow a 2D eonvolution model 
between the point spread function (PSF) and the tissue refleetivity funetion (TRF), i.e. the image to be 
recovered [|22ll . Speeifically, this results in r = Hx, where H G is a bloek cireulant with eirculant 

bloek (BCCB) matrix related to the 2D PSF of the system and x G represents the lexicographieally 
ordered tissue refleetivity funetion. We emphasize that this eonvolution model is based on the assumption 
of spatially invariant PSF. Although the PSF of an ultrasound system varies spatially and mostly in the 
axial direetion of the images, many settings of the ultrasound imaging system allow to attenuate this spatial 
variation in praetiee, sueh as the dynamie foealization of the reeeived eehoes, the multiple foeusing in 
emission and the time gain eompensation (TGC). For this reason, the PSF is eonsidered spatially invariant 
in our work, as it is the ease in part of the existing works on image deeonvolution in ultrasound imaging 
(e.g. unilloll). 

The objective of our paper is to propose a novel teehnique whieh is able to jointly aehieve US data 
volume reduetion and image quality improvement. In other words, the main idea is to combine the two 
frameworks of CS and deeonvolution applied to US imaging, resulting in the so-ealled eompressive deeon¬ 
volution (or CS deblurring) problem [l23l - [27ll . The combined direet model of joint CS and deconvolution 
is as follows: 


y = ^Hx -|- n (2) 

where the variables y, $, H, x and n have the same meaning as defined above. Inverting the model 
in Q will allow us to estimate the TRF x from the eompressed RF measurements y. 

To our knowledge, our work is the first attempt of addressing the eompressive deeonvolution problem 
in US imaging. In the general-purpose image proeessing literature, a few methods have been already 
proposed aiming at solving @ [|23l - l29ll . In [|25l l28l l29l . the authors assumed x was sparse in the direet 
or image domain and the PSF was unknown. In [|2^ [29ll . a study on the number of measurements lower 
bound is presented, together with an algorithm to estimate the PSF and x alternatively. The authors in [l25ll 
solved the eompressive deeonvolution problem using an £i-norm minimization algorithm by making use of 
the ”all-pole” model of the autoregressive proeess. In [|23ll24l, x was eonsidered sparse in a transformed 
domain and the PSF was supposed known. An algorithm based on Poisson singular integral and iterative 
curvelet thresholding was shown in ff^ . The authors in [|24l further combined the curvelet regularization 
with total variation to improve the performanee in [|23l. Finally, the methods in [|^ [27ll supposed the 
blurred signal r = Hx was sparse in a transformed domain and the PSF unknown. They proposed a 
eompressive deeonvolution framework that relies on a eonstrained optimization teehnique allowing to 
exploit existing CS reeonstruetion algorithms. 

In this paper, we propose a eompressive deeonvolution teehnique adapted to US imaging. Our solution 
is based on the alternating direetion method of multipliers (ADMM) [fSOllSTll and exploits two eonstraints. 
The first one, inspired from CS, imposes via an Zi-norm the sparsity of the RF image r in the transformed 
domain (Fourier domain, wavelet domain, etc). The seeond one imposes a priori information for the TRF 
X. Gaussian and Laplaeian statisties have been extensively explored in US imaging (see e.g. lITTl 1^ 1^ . 
Moreover, reeent results show that the Generalized Gaussian Distributed (GGD) is well adapted to model 
the TRF. Consequently, we employ herein the minimization of an /p-norm of x, covering all possible cases 
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ranging from 1 to 2 [|T9l 1201 [34ll . Similar to all existing frameworks, we eonsider the CS sampling matrix 
$ known. In US imaging, the PSF is unknown in practical applications. However, its estimation from the 
RF data as an initialization step for the deconvolution has been extensively explored in US imaging. In 
this paper, we adopted the approach in [f35]l in order to estimate the PSF further used to construct the 
matrix H. 

This paper is organized as follows. In section |I^ we formulate our problem as a convex optimization 
routine and propose an ADMM-based method to efficiently solve it. Supporting simulated and experimental 
results are provided in section showing the contribution of our approach compared to existing methods 


and its efficiency in recovering the TRF from compressed US data. The conclusions are drawn in IV 


II. Proposed Ultrasound Compressive Deconvolution Algorithm 
A. Optimization Problem Formulation 

1) Sequential approach: In order to estimate the TRF x from the compressed and blurred measurements 
y, an intuitive idea to invert the direct model in Q is to proceed through two sequential steps. The aim 
of the first step is to recover the blurred US RF image r = Hx from the compressed measurements y 
by solving the following optimization problem: 


min ||a|| i H- \\y — I 

2/i 


(3) 


where a is the sparse representation of the US RF image r in the transformed domain 'F, that is, 
r = Hx = Different basis have been shown to provide good results in the application of CS in US 
imaging, such as wavelets, waveatoms or 2D Fourier basis In this paper the wavelet transform has 
been employed. 

Once the blurred RF image, denoted by r, is recovered by solving the convex problem in Q, one can 
restore the TRF x by minimizing: 


min a 


1*11^ + 


Hx\ 


(4) 


The first term in Q aims at regularizing the TRF by a GGD statistical assumption, where p is related 
to the shape parameter of the GGD. In this paper, we focus on shape parameters ranging from 1 to 2 
(1 < p < 2), allowing us to generalize the existing works in US image deconvolution mainly based on 
Laplacian or Gaussian statistics ifTbl - fT^ . 

2) Proposed approach: While the sequential approach represents the most intuitive way to solve 
the compressive deconvolution problem, dividing a single problem into two separate subproblems will 
inevitably generate larger estimation errors as shown by the results in section III Therefore, we propose 
herein a method to solve the CS and deconvolution problem simultaneously. Similarly to [|^ . we formulate 
the reconstruction process into a constrained optimization problem explointing the relationship between 
(|3]) and @. 


min 

a:eIR^,aeK^ 


a 


ill +T 



1 

2p 


y - II 2 


(5) 


s.t. Hx = \Fa 

However, since our goal is to recover enhanced US imaging by estimation the TRF x, we further 
reformulate the problem above into a unconstrained optimization problem: 


min II ^ ^Hx 111 +a \\x\\l + II 2/ - ^Hx || I 

a;eR^ ^ 2/i 


( 6 ) 


The objective function in (|^ contains, in addition to the data fidelity term, two regularization terms. 
The first one aims at imposing the sparsity of the RF data Hx (i.e. minimizing the fi-norm of the target 
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image x convolved with a bandlimited function) in a transformed domain \I/. We should note that such an 
assumption has been extensively used in the application of CS in US imaging, see e.g. flU [6]-[8l |T3l |34| . 
Transformations such as 2D Fourier, wavelet or wave atoms have been shown to provide good results 
in US imaging. The second term aims at regularizing the TRF x and is related to the GGD statistical 
assumption of US images, see e.g. [fT9l l20l l34ll . 

We notice that our regularized reconstruction problem based on the objective function in @ is different 
from a typical CS reconstruction. Specifically, the objective function of a standard CS technique applied 
to our model would only contain the classical data fidelity term and an £i-norm penalty similar to the first 
term in @ but without the operator H. However, such a CS reconstruction is not adapted to compressive 
deconvolution, mainly because the requirements of CS theory such as the restricted isometry property 
might not be guaranteed [l26l . 

To solve the optimization problem in eq. (|^, we propose hereafter an algorithm based on the alternating 
direction method of multipliers (ADMM). 


B. Basics of Alternating Direction Method of Multipliers 

Before going into the details of our algorithm, we report in this paragraph the basics of ADMM. ADMM 
has been extensively studied in the areas of convex programming and variational inequalities, e.g., [f30l . 
The general optimization problem considered in ADMM framework is as follows: 


min f{u) + g{v) 

U^V 

s.t. Bu + Cv = b,u^ll,v^V 


(7) 


where U C W and V C M* are given convex sets, / : W —M and : V —?■ M are closed convex 
functions, B G and C G are given matrices and b G is a given vector. 

By attaching the Lagrangian multiplier A G to the linear constraint, the Augmented Lagrangian (AL) 
function of Q is 


C{u, V, A) = f{u) + g{v) — X^{Bu + Cv — b) 

+ ^ \\ Bu + Cv-b \\l ^ ^ 

where /? > 0 is the penalty parameter for the linear constraints to be satisfied. The standard ADMM 
framework follows the three steps iterative process: 

G argminC{u,v^, 

u&U 

G argminC{u^~^^,v,X^) (9) 

v&V 

yk+i = yfc _ + Cv^+^ - b) 

The main advantage of ADMM, in addition to the relative ease of implementation, is its ability to split 
awkward intersections and objectives to easy subproblems, resulting into iterations comparable to those 
of other first-order methods. 


C. Proposed ADMM parameterization for Ultrasound Compressive Deconvolution 

In this subsection, we propose an ADMM method for solving the ultrasound compressive deconvolution 
problem in Q. 

Using a trivial variable change, the minimization problem in Q can be rewritten as: 


min II w 111 +a ||a;||!? H- W V — Aa 

^ 2/i 


2 

2 


( 10 ) 
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where a = ^ ^Hx, w = a and A = Let us denote z = 


B = 


In 


,C' = 


0 


-H 


The reformulated problem in (10) 


ean fit the general ADMM framework in (|^ by ehoosing: /(a) = ^ 

^ and 6 = 0. Jjv G is the identity matrix. 


y-Aa III, g{z) =11 w ||i +a ||c^, 


The augmentea Lagrangian funetion of ([TO]) is given by 


C{a, z, A) = /(a) + g{z) — \^{Ba + Cz) 


+ — II Ba + Cz 
2 " 


( 11 ) 


where A G 


d27V 


stands for A = 


Ai 

A2 


, Ai G = 1,2). Aeeording to the standard ADMM iterative 


seheme, the minimizations with respeet to a and z will be performed alternatively, followed by the update 

of A. 


D. Implementation Details 

In this subsection, we detail each of the three steps of our ADMM-based compressive deconvolution 
method. While the following mathematical developments are given for the case when the regularization 
term for TRF x is equal to || x \\^ (adapted to US images), our approach using a generalized total variation 

be useful for other (medical) applications. 
w 


regularization is also detailed in Appendix A and may 
Step 1 consists in solving the z-problem, since z = 
subproblems. 


X 


, this problem can be further divided into two 


Step 1.1 aims at solving: 


w =argmin 


w 111 -(A^ 




W] 


+ — II ^ ~ II2 
2 " 


k • I, „ I, k-i K 

^ w =argmm || lu ||i +— || a — w - — 


k-l 




( 12 ) 


^ w — 




a 


k-l 


K 


k-l 


fl 


where prox stands for the proximal operator as proposed in Il36] - I38ll . The proximal operators of various 
kinds of functions including ||a;||p have been given explicitly in the literature (see e.g. ||39l). Basics about 
the proximal operator of ||a;||!’ are reminded in Appendix B. 


Step 1.2 consists in solving: 


x^ = argmin a\\ x ||^ —A 2 (f&a 


k-l 




Hx) 


+ - II 
2 " 


Hx III 


(13) 


For p equal to 2, the minimization in ([T3]) can be easily solved in the Fourier domain, as follows: 


X 


= [/3H^H + 2aJiv] ^ x 


(14) 


For 1 ^ p < 2, we propose to use the proximal operator to solve ( |13[ ). In this case, eq. 0 will be 
equivalent to 
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= argmin a || a; ||^ || 


.k-l 


Hx 


. A:—1 
'2 ||2 
2 




(15) 


Denoting h{x) = | || ^ — Hx — 


.k-l 

12 


/3 II 2 
^/c—1 


, we ean further approximate h{x) by 


h'{x^ ^)(a; — x^ H-|| a; — a; 

27 


i|2 


(16) 


where 7 > 0 is a parameter related to the Lipschitz eonstant [|40l and h'{x^~^) is the gradient of h{x) 
when X = x^~^. 

By plugging ( [T^ into ([T5]), we obtain: 


x^ ^argmin a\\x ||j^ +I3h'{x^ ^)(a; — x^ ^) 

tcSK^ 


+ 111 

27 


/S 


77 x^ ^argmin a || a; ||^ +— || a; — a; 
cceM^ 27 


k-l 


+ -ih'{x^-^) 


Aeeording to the definition of the proximal operator, we ean finally get 


(17) 


x^ ^ prox^^\\.\\v/^{x^ ^--ih\x^ ^)} (18) 

We should note that ( [T8] ) provides an approximate solution, thus resulting into an inexaet ADMM 
seheme. However, the convergenee of such inexact ADMM has been already established in [1^ 1411 l42l . 
Step 2 aims at solving: 


= argmin — \\ y - Aa -(A^-^)*(Sa + Cz’^) 

aeR^ 2/i 

+ ^ II Ba + Cz’^ ||2 

77 = i-A^A + (3In + 

y y 

+ (5w^ + (5^^Hx^) 


(19) 


The formula above is equivalent to solving xci N x N linear system or inverting an A^ x A^ matrix. 
However, since the sparse basis 'h considered is orthogonal (e.g. the wavelet transform), it can be reduced 
to solving a smaller M x M linear system or inverting an M x M matrix by exploiting the Sherman- 
Morrison-Woodbury inversion matrix lemma P3l : 


- ^A\/3Jm + ^2AAy^A ( 20 ) 

Pi Pi 

In this paper, without loss of generality, we considered the compressive sampling matrix <I) as a 
Structurally Random Matrix (SRM) ll4^ . Therefore, A was formed by randomly taking a subset of rows 
from orthonormal transform matrices, that is, A A* = Jm- As a consequence, there is no need to solve a 
linear system and the main computational cost consists into two matrix-vector multiplications per iteration. 
Step 3 consists in solving: 


Afc = + cz^) 

The proposed optimization routine is summarized in Algorithm 


( 21 ) 
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Algorithm 1 ADMM algorithm for Solving (6) 
Input: a°, A°, a, fi, (5 
1 : while not converged do 
2 : 

3 : ^ 

4 : <(— W^, 

5 : , x^, 

6: end while 
Output: X 


> update using 
> update x^ using ( [T4| ) or 


( 12 ) 


(18) 


> update using ( fT9 |) 

> update using 


III. Results 

The performanee of the proposed eompressive deeonvolution method are evaluated on several simulated 
and experimental data sets. First, we test our algorithm on a modified Shepp-Logan phantom eontaining 
speckle noise to eonfirm that the /p-norm regularization term is more adapted to US images than the 
generalized TV used in [l26l . The approach in [l26l is referred as CD Amizie hereafter. Second, we 
give the results of our algorithm for different /p-norm optimizations on simulated US images, showing 
the superiority of our approach over the intuitive sequential method explained in seetion |I^ Finally, 
eompressive deeonvolution results on two in vivo ultrasound images are presented. Moreover, a eomparison 
between our approaeh and CD Amizie on the standard Shepp-Logan phantom is provided in Appendix 
C. 


A. Results on modified Shepp-Logan phantom 

We modified the Shepp-Logan phantom in order to simulate the speekle noise that degrades in praetiee 
the US images. For this, we followed the procedure classically used in US imaging [|45l. First, seatterers 
at uniformly random loeations have been generated, with amplitudes distributed aeeording to a zero-mean 
generalized Gaussian distribution (GGD) with the shape parameter set to 1.3 and the seale parameter equal 
to 1. The scatterer amplitudes were further multiplied by the values of the original Shepp-Logan phantom 
pixels loeated at the elosest positions to the seatterers. The resulting image, mimieking the tissue refleetivty 
funetion (TRF) in US imaging, is shown in Figj^a). The blurred image in Figj^b) was obtained by 2D 
eonvolution between the TRF and a spatially invariant PSF generated with Field II [146]|, a state-of-the-art 
simulator in US imaging. It eorresponds to a 3.5 MHz linear probe, sampled in the axial direetion at 20 
MHz. The eompressive measurements were obtained by projeeting the blurred image onto SRM and by 
adding a Gaussian noise eorresponding to a SNR of 40 dB. 

The results were quantitatively evaluated using the standard peak signal-to-noise ratio (PSNR) and the 
Struetural Similarity (SSIM) [|47l . PSNR is defined as 

NL‘^ 

PSNR = lOlogio-r -^ (22) 

II a: — a: 11 

where x and x are the original and reconstrueted images, and the eonstant L represents the maximum 
intensity value in x. SSIM is extensively used in US imaging and defined as 


SSIM 


(2,Pxdx T Cl)(2(7 XX + C 2 ) 

i.dl +l4 + Ci){al + al + C2) 


(23) 


where x and x are the original and reeonstrueted images, dx, ^x and are the mean and varianee 
values of x and x, a^x is the eovarianee between x and x; ci = {kiLfi and C 2 = (/c 2 -^)^ are two variables 
aiming at stabilizing the division with weak denominator, L is the dynamie range of the pixel-values and 
ki, k 2 are eonstants. Herein, L = 1, ki = 0.01 and ^2 = 0.03. 














(e) (f) (g) 


Fig. 1: Reconstruction results for SNR = 40dB and a CS ratio of 0.6. (a) Modified Shepp-Logan phantom 
containing random scatterers (TRF), (b) Degraded image by convolution with a simulated US PSF, 
(c) Reconstruction result with CD Amizic, (d) Reconstruction result with the proposed method using 
a generalized TV prior (ADMM GTV), (e, f, g) Reconstruction results with the proposed method using 
an /p-norm prior, for p equal to 1.5, 1.3 and 1 (ADMM_L1.5, ADMM_L1.3 and ADMM_L1). 


Reconstruction results for a CS ratio of 0.6 are shown in Figj^ They were obtained with: the recent 
compressive deconvolution technique reported in ff2^ (referred as CD_Amizic), the proposed method using 
the generalized TV prior (denoted by ADMM GTV) and the proposed method using the /p-norm for p 
equal to 1.5, 1.3 and 1 (denoted respectively by ADMM_L1.5, ADMM_L1.3 and ADMM Ll). All the 
hyperparameters were set to their best possible values by cross-validation. For CD Amizic, {/3, a, p, r} = 
{10^, 1,10^, 10^}. For ADMM_GTV {p, a, (5} = 2 x 10“^, 10^} and for the proposed method with 

/p-norms, {/r, a,/9,7} = {10“^,2 x 10“^, 10^,3 x 10“^} . The quantitative results for different CS ratios 
are regrouped in Tablej^ They confirm that the Zp-norm is better adapted to recover the TRF in US imaging 
than the generalized TV. The difference between the two priors is further accentuated when the CS ratio 
decreases. 

Keeping in mind that the generalized TV prior is not well suited to model the TRF in US imaging, we 
did not use CD Amizic in the following sections dealing with simulated and experimental US images. 
Moreover, the proposed method was only evaluated in its Zp-norm minimization form. 

TABLE I: Quatitative results for the modified Shepp-Logan phantom with US speckle (SNR = 40dB) 


CS ratios 


CDAmizic 

ADMM GTV 

ADMM L1.5 

ADMM L1.3 

ADMM L1 

80% 

PSNR 

SSIM 

30.82 

83.24 

31.II 

85.03 

32.23 

86.44 

32.32 

88.77 

32.05 

87.70 

60% 

PSNR 

SSIM 

29.68 

74.58 

29.83 

77.83 

31.27 

82.26 

31.50 

86.03 

31.32 

85.64 

40% 

PSNR 

SSIM 

26.76 

43.43 

28.11 

61.46 

29.58 

73.88 

30.04 

79.95 

30.12 

81.75 

20% 

PSNR 

SSIM 

20.22 

8.35 

21.53 

12.77 

26.81 

51.70 

27.29 

62.93 

28.20 

72.34 


B. Results on simulated ultrasound images 

In this section, we compared the compressive deconvolution results with our method to those obtained 
with a sequential approach. The latter recovers in a first step the blurred US image from the CS measure- 
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(d) (e) (f) 

Fig. 2: Simulated US image and its eompressive deeonvolution results for a CS ratio of 0.4 and a SNR of 
40 dB. (a) Original tissue refleetivity funetion, (b) Simulated US image, (e) Results using the sequential 
method, (d, e, f) Results with the proposed method for p equal to 1.5, 1.3 and 1 respeetively. 


ments and reeonstruets in a seeond step the TRF by deeonvolution. 

Two ultrasound data sets were generated, as shown in Figures and They were obtained by 2D 
eonvolution between spatially invariant PSFs and the TRFs. For the first simulated image, the same PSF 
as in the previous seetion was simulated and the TRF eorresponds to a simple medium representing a 
round hypoeehoie inelusion into a homogeneous medium. The seatterer amplitudes were random variables 
distributed aeeording to a GGD with the shape parameter set to 1. The seeond data set is one of the 
examples proposed by the Field II simulator [l46ll . mimieking a kidney tissue. The PSF was also generated 
with Field II eorresponding to a 4 MHz eentral frequeney and an axial sampling frequeney of 40 MHz. It 
eorresponds to a foealized emission (the PSF was measured at the foeal point) with a simulated linear probe 
eontaining 128 elements. The shape parameter of the GGD used to generate the seatterer amplitudes was 
set to 1.5 and the number of seatterers was eonsidered suffieiently large (10®) to ensure fully developped 
speekle. In both experiments, the eompressed measurements were obtained by projeeting the RF images 
on SRM, aiming at redueing the amount of data available. 

With the sequential approaeh, YALLl [l42ll was used to proeess the CS reeonstruetion following the 
minimization in eq. Q. The deeonvolution step was proeessed using the Forward-Baekward Splitting 
method [|48l |49l . Both the CS reeonstruetion and the deeonvolution proeedures were performed with the 
same priors as the proposed eompressive deeonvolution approaeh. 

The algorithm stops when the eonvergenee eriterion || || / || x^~^ ||< le“^ is satisfied. In order 

to highlight the influenee of these hyperparameters on the reeonstruetion results, we eonsider the simulated 
US image in Fig. The PSNR values obtained while varying the values of these hyperparameters are 
shown in Fig. From Fig. Q one ean observe that the best results are obtained for small values of p, 
eorresponding to an important weight given to the data attaehment term. The best value of a is the one 
providing the best eompromise between the two prior terms eonsidered in eq. promoting minimal 
fi-norm of Hx in the wavelet domain and GGD statisties for x. The ehoiee of j3 and 7 parameters, used 
in the augmented Lagrangian funetion and in the approximation of the £p-norm proximal operator, have 
an important impaet on the algorithm eonvergenee. Moreover, one may observe that for a given range of 
values, the ehoiee of 7 has less impaet on the quality of the results than the other three hyperparameters. 
Despite different optimal values for eaeh CS ratio, in the results presented through the paper, we eonsidered 
their values fixed for all the CS ratios. The hyperparameters with our approaeh were set to {/r, a, /), 7} = 
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(d) (e) (f) 

Fig. 3: Simulated kidney image and its eompressive deeonvolution results for a CS ratio of 0.2 and a SNR 
of 40dB. (a) Original tissue refleetivity funetion, (b) Simulated US image, (e) Results using the sequential 
method, (d, e, f) Results with the proposed method for p equal to 1.5, 1.3 and 1 respeetively. 


—o— CSra1io=0.8 —®— CSratio=0.6 — ®— CSratio=0.4 — ^CSratio=0.2| 



■-7 -6 -5 -4 -3 -2 -3 - 2.5 -2 - 1.5 -1 - 0.5 0 0.5 1 





Fig. 4: The impaet of hyperparameters on the performanee of proposed algorithm on Figure. 















11 


TABLE II : Quantitative results for simulated US images (SNR = 40dB) 


CS 

Ratios 


Sequential 

Proposed 

(fi.s) 

Proposed 

ih.s) 

Proposed 

(fi) 

Figure 

) 

80% 

PSNR 

SSIM 

26.50 

75.01 

24 

73 

.74 

.91 

25.29 

77.66 

26.82 

79.45 

60% 

PSNR 

SSIM 

25.96 

68.59 

24.44 

69.37 

24.74 

74.72 

26.03 

76.26 

40% 

PSNR 

SSIM 

23.38 

47.60 

24.21 

62.58 

24.57 

71.86 

25.28 

72.78 

20% 

PSNR 

SSIM 

21.10 

36.07 

23.72 

50.34 

24.42 

66.48 

24.77 

70.44 

Figure 


80% 

PSNR 

SSIM 

26.06 

45.99 

25 

56 

.81 

26.72 

56.84 

26.69 

56.71 

60% 

PSNR 

SSIM 

25.44 

38.86 

26.38 

54.14 

26.31 

53.90 

26.29 

53.80 

40% 

PSNR 

SSIM 

25.37 

34.61 

25.89 

50.22 

25.95 

50.51 

25,97 

50.61 

20% 

PSNR 

SSIM 

24.96 

30.89 

25.22 

41.41 

25.20 

41.32 

25.12 

40.97 


2 X 10“^, 1,10“^} for the round cyst image and {/i, a, /3, 7} = 2 x 10“^, 1 x 10^, 10“'^} for 

the simulated kidney image. 

The quantitative results in Table show that the proposed method outperforms the sequential approach, 
for all the CS ratios and values of p considered. They confirm the visual impression given by Figures 
and 1^ We should remark that for the first simulated data set, the /i-norm gives the best result. This may 
be explained by the simple geometry of the simulated TRF, namely its sparse appearance. The second 
data set, more realistic and more representative of experimental situations, shows the interest of using 
different values of p. It confirms the generality interest of the proposed method, namely its flexibility in 
the choice of TRF priors. 


C. In vivo study 

In this section, we tested our method with two in vivo data sets. The experimental data were acquired 
with a 20 MHz single-element US probe on a mouse bladder (first example) and kidney (second example). 
Unlike the simulated cases studied previously, the PSF is not known in these experiments and has to be 
estimated from the data. In this paper, the PSF estimation method presented in [[^ has been adopted. 
The PSF estimation adopted is not iterative and the computational time for this pre-processing step is 
negligible compared to the reconstruction process. The compressive deconvolution results are shown in 
Figures and for different CS ratios. 

Given that the true TRF is not known in experimental conditions, the quality of the reconstruction 
results is evaluated using the contrast-to-noise ratio (CNR) [l50ll . defined as 


CNR = 


\dl ~ d2\ 


(24) 


where and p 2 are the mean of pixels located in two regions extracted from the image while cxi and a 2 
are the standard deviations of the same blocks. The two regions selected for the computation of the CNR 
are highlighted by the two red rectangles in Figures [^a) and |^a). Table. Ill gives the CNR assessment 
for these two in vivo data sets with different CS ratios and p values. Given the sparse appearance of the 
bladder image in Fig. |^a), the best result was obtained for p equal to 1. However, the complexity of the 
tissue structures in the kidney image in Fig. [^results into better results for p larger than 1. Nevertheless, 
both the visual impression and the CNR results show the ability of our method to both recover the image 
from compressive measurements and to improve its contrast compared to the standard US image. In 
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Fig. 5: From left to right, the original in vivo image and its eompressive deeonvolution results for CS 
ratios of 1, 0.8, 0.6 and 0.4 respeetively with p = 1. 



(a) (b) (c) (d) (e) 


Fig. 6: From left to right, the original in vivo image and its eompressive deeonvolution results for CS 
ratios of 1, 0.8, 0.6 and 0.4 respeetively with p = 1.5. 


partieular, we may remark the improved eontrast of the struetures inside the kidney on our reeonstrueted 
images eompared to the original one. 


TABLE III: Cb 

IR assess 

ment for in vivo data 

Figure 

Original CNR 

p values 

CS ratios 

100% 80% 60% 40% 

Fig.6 

1.106 

p= 1 

p = 1.5 

1.748 1.546 1.367 1.333 

1.690 1.424 1.304 1.287 

Fig.7 

1.316 

p= 1 

p = 1.5 

2.373 2.162 1.895 1.434 
2.317 2.082 1.905 1.451 


IV. Conclusion 

This paper introdueed an ADMM-based eompressive deeonvolution framework for ultrasound imaging 
systems. The main benefit of our approaeh is its ability to reeonstruet enhaneed ultrasound RF images from 
eompressed measurements, by inverting a linear model eombining random projeetions and 2D eonvolution. 
Compared to a standard eompressive sampling reeonstruetion that operates in the sparse domain, our 
method solves a regularized inverse problem eombining the data attaehment and two prior terms. One 
of the regularizers promotes minimal £i-norm of the target image transformed by 2D eonvolution with 
a bandlimited ultrasound PSF. The seeond one is seeking for imposing GGD statisties on the tissue 
refleetivity funetion to be reeonstrueted. Simulation results in Appendix C on the standard Shepp-Logan 
phantom show the superiority of our method, both in aeeuraey and in eomputational time, over a reeently 
published eompressive deeonvolution approaeh. Moreover, we show that the proposed joint CS and 
deeonvolution approaeh is more robust than an intuitive teehnique eonsisting of first reeonstrueting the RF 
data and seeond deeonvolving it. Finally, promising results on in vivo data demonstrate the effeetiveness of 
our approaeh in praetieal situations. We emphasize that the 2D eonvolution model may not be valid over 
the entire image beeause of the spatially variant PSF. While in this paper we foeused on eompressive image 
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deconvolution based on spatially invariant PSF, a more complicated global model combining several local 
shift invariant PSFs represents an interesting perspective of our approach. Our future work includes: I) the 
consideration of blind deconvolution techniques able to estimate (update) the spatially variant or invariant 
PSF during the reconstruction process, II) automatic techniques for choosing the optimal parameter p used 
to regularize the tissue reflectivity function, III) extend our method to 3D ultrasound imaging, IV) evaluate 
other existing setups to obtain the random compressed measurements further adapted to accelerate the 
frame rate instead of only reducing the amount of acquired data, V) consider the compressed deconvolution 
of temporal image sequences by taking advantage of the information redundancy and by including in our 
model the PSF frame-to-frame variation caused by strong clutters in in vivo scenarios, VI) evaluate our 
approach on more experimental data. 


Appendix A 

Proposed compressive deconvoeution with generalized TV 

Although the generalized total variation (TV) used in is not suitable for ultrasound images, it may 
have great interest in other (medical) application dealing with piecewise constant images. As suggested 
in [|26l , the generalized TV is given by: 


( 25 ) 

deD i 

where o{d) G {1,2} denotes the order of the difference operator Af(a;), 0 < p < 1, and d e D = 
{h,v,hh,vv,hv}. A^{x) and A}(a;) correspond, respectively, to the horizontal and vertical first order 
differences, at pixel i, that is, A^{x) = Ui — ui^i) and A}(a;) = Ui — Ua{i), where l{i) and a{i) denote 
the nearest neighbors of i, to the left and above, respectively. The operators Af^(a;), A™(a;), A^'^{x) 
correspond, respectively, to horizontal, vertical and horizontal-vertical second order differences, at pixel 
i. 

Replacing the £p-norm by the generalized TV in our compressive deconvolution scheme results in a 
modified x update step, that turns in solving: 


X =argrmn a 


E 2'^°“ El AW 




- - Hx) + ^ II - Hx 


Similarly to the first step of the method in II261I . the equation above can be solved iteratively by: 


x^’^ = 


-1 


d 


(26) 


where I is the iteration number in the process of updating x, is a diagonal matrix with entries A*^ 
is the convolution matrix (BCCB matrix) of the difference operator Af(-) and which 

is updated iteratively by: 


V 


d,i 




2 


When a stopping criterion is met, we can finally get an update of x. 


( 27 ) 
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Appendix B 
Proximal Operator 

The proximal operator of a function / is defined for e by: 

proxf{x^) = argmin f{x) + - || t II 2 (28) 

£ceM^ 2 

When f = K\ x\^, the corresponding porximal operator has been given by [[391: 

proxK\x\p{x^) = sign{x^)q (29) 

where q ^ 0 and 

q + pKq^~^ = |a;°| (30) 

It is obvious that the proximal operator of K |a;| is a soft thresholding, which is equal to: 

proxK\x\{,x^) = max [\x^\- (31) 

When p 7 ^ 1, we used Newton’s method to obtain its numerical solution, i.e. the value of q. 

Appendix C 

Comparison with the compressive deconvolution method in [|2^ 

In this appendix we show an experiment aiming to evaluate the performance of the proposed approach 
compared to the one in [|2i6l , denoted by CD_Amizic herein. The comparison results are obtained on 
the standard 256 x 256 Shepp-Logan phantom. The measurements have been generated in a similar 
manner as in [j26l, i.e. the original image was normalized, degraded by a 2D Gaussian PSF with a 5- 
pixel variance, projected onto a structured random matrix (SRM) to generate the CS measurements and 
corrupted by an additive Gaussian noise. We should remark that in [|26l the compressed measurements 
were originally generated using a Gaussian random matrix. However, we have found that the reconstruction 
results with CD Amizic are slightly better when using a SRM compared to the PSNR results reported 
in [26]. Both methods were based on the generalized TV to model the image to be estimated and the 
3-level Haar wavelet transform as sparsifying basis 'P. With our method, the hyperparameters were set to 
/)} = 10“^, 10^}. The same hyperparameters as reported in [|26ll were used for CD_Amizic. 

Both algorithms based on the non-blind deconvolution (PSF is supposed to be known) and used the same 
stopping criteria. 

FigjT] shows the original Shepp-Logan image, its blurred version and a series of compressive decon¬ 
volution reconstructions using both our method and CD Amizic for CS ratios running from 0.4 to 0.8 
and a SNR of 40 dB. Table jlV] regroups the PSNRs obtained with our method and with CD Amizic 
for two SNRs and for four CS ratios from 0.2 to 0.4. In each case, the reported PSNRs are the mean 
values of 10 experiments. We may observe that our method outperforms CD Amizic in all the cases, 
allowing a PSNR improvement in the range of 0.5 to 2 dB. Moreover, Fig|^ shows the computational 
times with CD Amizic and the proposed method, obtained with Matlab implementations (for CD Amizic, 
the original code provided by the authors of ll26ll has been employed on a standard desktop computer 
(Intel Xeon CPU E5620 @ 2.40GHz, 4.00G RAM). We notice that our approach is less time consuming 
than CD_Amizic for all the CS ratios considered. 
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(a) (b) (c) (d) 



(e) (f) (g) (h) 


Fig. 7: Shepp-logan image and its compressive deconvolution results for a SNR of 40dB. (a) Original 
image, (e) Blurred image, (b,c,d) Compressive deconvolution results with CD_Amizic for CS ratios of 
0.8, 0.6 and 0.4, (f,g,h) Compressive deconvolution results with the proposed method for CS ratios of 
0.8, 0.6 and 0.4. 


400 


350 


300 
^ 250 

P 

D) 200 

'c 

c 150 

cr 

100 

50 

0 


Fig. 8: Mean reconstruction running time for 10 experiments conducted for each CS ratio for a SNR of 
40 dB. 
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TAB LE IV: PSNR assessment for Shepp-Logan phan tom 


SNR 

CS ratios 

20% 

40% 

60% 

80% 

40dB 

CDAmizic 
Proposed method 

23.04 

24.09 

24.88 

25.38 

25.30 

26.26 

25.51 

26.91 

30dB 

CD_Amizic 
Proposed method 

22.61 

23.92 

24.05 

25.12 

24.40 

25.82 

24.55 

26.33 
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